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Enzymes utilize protein architectures to create highly specialized structural motifs that can greatly enhance 
the rates of complex chemical transformations. Here we use experiments, combined with ab initio simulations 
that exactly include nuclear quantum effects, to show that a triad of strongly hydrogen bonded tyrosine 
residues within the active site of the enzyme ketosteroid isomerase (KSI) facilitates quantum proton delo¬ 
calization. This delocalization dramatically stabilizes the deprotonation of an active site tyrosine residue, 
resulting in a very large isotope effect on its acidity. When an intermediate analog is docked, it is incorpo¬ 
rated into the hydrogen bond network, giving rise to extended quantum proton delocalization in the active 
site. These results shed light on the role of nuclear quantum effects in the hydrogen bond network that sta¬ 
bilizes the reactive intermediate of KSI, and the behavior of protons in biological systems containing strong 
hydrogen bonds. 


I. INTRODUCTION 

While many biological processes can be well described 
with classical mechanics, there has been much interest 
and debate as to the role of quantum effects in biological 
systems ranging from photosynthetic energy transfer, to 
photoinduced isomerization in the vision cycle and avian 
magnetoreception 1 . For example, nuclear quantum ef¬ 
fects, such as tunneling and zero-point energy (ZPE), 
have been observed to lead to kinetic isotope effects of 
greater than 100 in biological proton and proton-coupled 
electron transfer processes 2,3 . However, the role of nu¬ 
clear quantum effects in determining the ground state 
thermodynamic properties of biological systems, mani¬ 
festing as equilibrium isotope effects, has gained signifi¬ 
cantly less attention 4 . 

Ketosteroid isomerase (KSI) possesses one of the high¬ 
est enzyme unimolecular rate constants and is thus 
considered a paradigm of proton transfer catalysis in 
enzymology 6 12 . KSI’s remarkable rate is intimately con¬ 
nected to the formation of a hydrogen bond network in 
its active site (Fig. la) which acts to stabilize a charged 
dienolate intermediate, lowering its free energy by ~11 
kcal/mol relative to solution (Fig. SI) 7 . This extended 
hydrogen bond network in the active site links the sub¬ 
strate to Aspl03 and Tyrl6, with the latter further hy¬ 
drogen bonded to Tyr57 and Tyr32 as shown in Fig. la. 

The mutant KSF 040 ^ preserves the structure of the 
wild-type enzyme while mimicking the protonation state 
of residue 40 in the intermediate complex (Fig. lb), 
therefore permitting experimental investigation of an 
intermediate-like state of the enzyme 7,13 15 . Experi¬ 
ments have identified that, in the absence of an in¬ 
hibitor, one of the residues in the active site of KSI D40Ar 
is deprotonated 13 . Although one might expect the car- 
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FIG. 1. KSI-intermediate and KSI 1540 ^-inhibitor complex. 
Schematic depiction of (a) the KSI-intermediate complex dur¬ 
ing the catalytic cycle (Fig. SI) and (b) a complex between 
Kgi^ov an q pp eno ^ an inhibitor which acts as an interme¬ 
diate analog. Both the intermediate and inhibitor are stabi¬ 
lized by a hydrogen bond network in the active site of KSI. 
(c) Image of KSI D40iV with the tyrosine triad enlarged and 
the atoms 016, H16, 032, H32 and 057 labeled (shown with 
Tyr57 deprotonated) 5 . 


boxylic acid of Asp 103 to be deprotonated, the combi¬ 
nation of recent 13 C NMR and UV-Vis experiments has 
shown that the ionization resides primarily on the hy¬ 
droxyl group of Tyr57, which possesses an anomalously 
low pK a of 6.3 zb 0.1 13 . Such a large tyrosine acidity 
is often associated with specific stabilizing electrostatic 
interactions (such as a metal ion or cationic residue in 
close proximity), which is not the case here suggesting 
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an additional stabilization mechanism is at play 16 . 

One possible explanation is suggested by the close 
proximity of the oxygen atoms (O) on the side chains 
of the adjacent residues Tyrl6 (016) and Tyr32 (032) 
to the deprotonated O on Tyr57 (057, Fig. lc) 5 . In sev¬ 
eral high-resolution crystal structures, these distances are 
found to be around 2.6 A 5 ’ 15,17 , which are much shorter 
than those observed in hydrogen bonded liquids such as 
water, where 0-0 distances are typically around 2.85 A. 
Such short heavy atom distances are only slightly larger 
than those typically associated with low-barrier hydro¬ 
gen bonds (LBHB) 18-20 , where extensive proton sharing 
is expected to occur between the atoms. In addition, at 
these short distances the proton’s position uncertainty 
(de Broglie wavelength) becomes comparable to the O- 
O distance, indicating that nuclear quantum effects could 
play an important role in stabilizing the deprotonated 
residue (Fig lc). In this work, we demonstrate how nu¬ 
clear quantum effects determine the properties of pro¬ 
tons in KSI D40N ’s active site hydrogen bond network in 
the absence and presence of an intermediate analog by 
combining ab initio path integral simulations and isotope 
effect experiments. 
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FIG. 2. Experimental UV-Vis titration curves of KSI^ 40 ^ 
(circles) and tyrosine (triangles) in buffered H 2 O (red) and 
D 2 O (blue). Fractional ionization of the phenolic moiety is 
monitored by measuring the change in absorption at 300 nm. 
All data sets were well fit to a sigmoid function, admitting 
four pK a values and two A pK a values, the acidity isotope 
effect. 


II. RESULTS AND DISCUSSIONS 

A. Isotope substitution experiments reveal large isotope 
effect on acidity 

To assess the impact of nuclear quantum effects on the 
anomalous acidity of Tyr57, we measured the isotope ef¬ 
fect on the acid dissociation constant upon substituting 
hydrogens (H) in the hydrogen bond network with deu¬ 
terium (D). Because tyrosinate absorbs light at 300 nm 
more intensely than tyrosine, titration curves were gen¬ 
erated by recording UV spectra of KSI D40Ar at different 
pU s (where L is H or D) 16 . These experiments (Fig. 2 
and Table SI) reveal a change in pK a upon H/D substitu¬ 
tion (A pK^ SI ) of 1.1 zb 0.14 for tyrosine in the KSI 1140 ^ 
active site. This pK a isotope effect is much larger than 
that observed for tyrosine in solution (A pK^° l = 0.53 
d= 0.08), and is also, to the best of our knowledge, the 
largest recorded pK a isotope effect 21 . Changes in static 
equilibrium properties, such as the pK a , upon isotope 
substitution arise entirely from the quantum mechanical 
nature of nuclei. Such a large excess isotope effect, de¬ 
fined as A A pK a = A pK^ SI — A pK^° l , of 0.57 zb 0.16 
thus indicates that the tyrosine triad in the active site of 
KSi^ 40Ar exhibits much larger nuclear quantum effects 
than those observed for tyrosine in aqueous solution. 

Given the possible role of nuclear quantum effects in 
KSI D40Ar , can one estimate how much the quantum na¬ 
ture of protons changes the acidity of Tyr57 compared 
to a situation in which all the nuclei in the enzyme ac¬ 
tive site were classical (m 00 )? In the quasi-harmonic 
limit one can show that the pK a varies as the inverse 
square root of the particle mass (m) 22,23 . By using this 


relation and the experimental pL values, we can extrap¬ 
olate to the classical (m -A 00 ) limit which yields that 
the pK a of Tyr57 in KSE 040 ^ would be 10.1 zb 0.5 if the 
hydrogens were classical particles (Fig. S2). Relative to 
the observed pK a of 6.3 zb 0.1, this implies nuclear quan¬ 
tum effects lower the pK a or Tyr57 by 3.8 zb 0.5 units: 
an almost four orders of magnitude change in the acid 
dissociation constant. 


B. Ab initio path integral simulations of KSI^ 40 ^ 

To provide insights into the molecular origins of the nu¬ 
clear and electronic quantum effects which stabilize the 
deprotonated Tyr57 residue, we performed simulations of 
KSI D40N . To treat the electronic structure in the active 
site we performed ab initio molecular dynamics (AIMD) 
simulations using a QM/MM approach 24 27 in which the 
QM region was treated by density functional theory at 
the B3LYP-D3 level 28,29 (see Methods). These simula¬ 
tions allow for bond breaking and formation as dictated 
by the instantaneous electronic structure rather than pre¬ 
defined bonding rules. 

AIMD simulations are typically performed treating the 
nuclei as classical particles. However, a classical treat¬ 
ment of the nuclei would predict that the pK a would not 
change upon isotope substitution. Nuclear quantum ef¬ 
fects can be exactly included in the static equilibrium 
properties for a given description of the electronic struc¬ 
ture using the path integral formalism of quantum me¬ 
chanics, which exploits the exact mapping of a system of 
quantum mechanical particles onto a classical system of 
ring polymers 30 33 . We combined this formalism with on- 
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the-fly electronic structure calculations and performed ab 
initio path integral molecular dynamics (AI-PIMD) sim¬ 
ulations of KSI^ 40 ^. These simulations treat both the 
nuclear and electronic degrees of freedom quantum me¬ 
chanically in the active site QM region, and also incor¬ 
porate the fluctuations of the protein and solvent envi¬ 
ronment in the MM region. The simulations consisted 
of between 47 and 68 QM atoms and more than 52,000 
MM atoms describing the rest of the protein and sol¬ 
vent (Table S2). These simulations, which until recently 
would have been computationally prohibitive, were made 
possible by accelerating the PIMD convergence using a 
generalized Langevin equation 34 , utilizing new methods 
to accelerate the extraction of isotope effects 22 and ex¬ 
ploiting graphical processing units (GPUs) to perform 
efficient electronic structure theory evaluations via an 
interface to the TeraChem code 35,36 . Such a combina¬ 
tion yielded almost three orders of magnitude speed-up 
compared to existing AI-PIMD approaches, allowing 1.1 
ps/day of simulation to be obtained using 6 NVIDIA Ti¬ 
tan GPUs. We have recently shown that AI-PIMD sim¬ 
ulations using the B3LYP-D3 functional give excellent 
predictions of isotope effects in water, validating such a 
combination for the simulation of isotope effects in hy¬ 
drogen bonded systems 37 . 


C. Quantum delocalization of protons in KSI^ 40 ^ 

The excess isotope effect, A A pK a , obtained from our 
simulations (see SI Methods Ic) was 0.50 db 0.03, in ex¬ 
cellent agreement with the experimental value of 0.57 d= 
0.16. The average distance between 057 and the ad¬ 
jacent 016 and 032 atoms obtained in our simulations 
were 2.56 A and 2.57 A, with a standard deviation in 
both cases of 0.09 A. The distribution of distances be¬ 
tween 016 and 057 explored in the simulation is shown 
in the Fig. 4a inset. These average 0-0 distances are 
slightly smaller than (and within the margin of error of) 
those in the starting crystal structure (~ 2.6 A) 5 . As we 
will discuss below, the close proximity of the neighboring 
016 and 032 groups plays a crucial role in the origins of 
the observed isotope effect. 

Figure 3a-c shows snapshots from AIMD simulations 
in which the nuclei are treated classically (Fig. 3a) or 
quantum mechanically using the path integral formalism 
(AI-PIMD, Fig. 3b and c), while videos of the simu¬ 
lation trajectories are provided in SI Videos 1-3. For 
the quantum simulations the H16 and H32 protons are 
shown as their full ring polymers, which arise from the 
path integral quantum mechanics formalism. The spread 
of the ring polymer representing each proton is related to 
its de Broglie wavelength (quantum mechanical position 
uncertainty) 38,39 . The uncertainty principle dictates that 
localization of a quantum mechanical particle increases 
its quantum kinetic energy. The protons will thus at¬ 
tempt to delocalize, i.e. spread their ring polymers, to 
reduce this energetic penalty. The resulting proton po¬ 


sitions in Fig. 3b arise from the interplay between the 
chemical environment, such as the covalent O-H bond, 
which acts to localize the proton and the quantum kinetic 
energy penalty that must be paid to confine a quantum 
particle. Inclusion of nuclear quantum effects thus allows 
the protons to delocalize between the hydroxyl oxygens to 
mitigate the quantum kinetic energy penalty (Fig. 3b), 
which is not observed classically (Fig. 3a). Confinement 
of D, which due to its larger mass has a smaller position 
uncertainty, leads to a much less severe quantum kinetic 
energy penalty and hence less delocalization (Fig. 3c). 

To characterize the degree of proton delocalization we 
define a proton sharing coordinate vx = dox,HX ~ 
do57,HX, where dox,HX is the distance of proton HX 
from oxygen atom OX and X=16 or 32. Hence ux = 0 
corresponds to a proton that is equidistant between the 
oxygen atoms of TyrX and Tyr57, while a positive value 
indicates proton transfer to Tyr57 from TyrX. Figures 
3 d-f show the probability distribution along the proton 
sharing coordinates v\§ and z / 32 for classical nuclei and 
quantum nuclei for H and D, respectively. The free ener¬ 
gies along zq 6 and z / 32 are provided in Fig. S3. In the clas¬ 
sical AIMD simulation, H16 and H32 remain bound to 
their respective oxygens throughout the simulation (zqg 
and z /32 are negative) with Tyr57 ionized 99.96% of the 
time (Fig. 3d). However, upon including nuclear quan¬ 
tum effects (AI-PIMD simulations) there is a dramatic 
increase in the range of values v\§ and z /32 can explore 
(Fig. 3e). In particular, the probability that Tyr57 is 
protonated (yx > 0) increases by about 150-fold for H 
upon including quantum effects (Fig. 3e), with the pro¬ 
ton “hole” equally shifted onto the adjacent Tyrl 6 or 
Tyr32 residues. Proton transfers between the residues 
are observed frequently (SI Videos 2 and 3) with site 
lifetimes on the order of 60 and 200 fs in the H and D 
simulations, respectively. Although PIMD simulations 
exactly include nuclear quantum effects for calculating 
static properties, they do not allow rigorous extraction 
of time-dependent properties; nevertheless, they offer a 
crude way to assess the time-scale of the proton motion. 
The frequent transfers observed are also consistent with 
Fig. 3e and f, which show a monotonic decrease in the 
probability along both zqg and z/ 32 , i.e. although the pro¬ 
ton transferred state is lower in probability, the proton 
transfer process along each of the proton sharing coor¬ 
dinates contains no free energy barrier (Fig. S3) and is 
thus kinetically fast. 

As an experimental counterpart, we used chemical 
shifts of 13 C^-Tyr labeled KSF 040 ^ as a measure of frac¬ 
tional ionization of each Tyr residue (SI Methods lb) 15 . 
This analysis yielded values of 79% for the Tyr57 ioniza¬ 
tion for H and 86 % for D compared with simulated values 
of 94.2% and 98.3% (± 0.3%), respectively. This repre¬ 
sents good quantitative agreement, since the population 
difference amounts to a difference in the relative free en¬ 
ergy between experiment and theory of 0.7 kcal/mol - 
an error which is within the expected accuracy of the 
electronic structure approach employed. In addition, the 
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FIG. 3. Delocalized protons in the active site of KSI D40iV from AIMD and AI-PIMD simulations. Snapshots of the active site 
of KSI D40iV (top panels) and probability distribution along the proton sharing coordinates viq and ^32 (bottom panels) when 
the nuclei are treated classically (Classical) or quantum mechanically for H (QM H) and D (QM D). In the top panels, green, 
red and white represent C, O and H atoms. The blue-grey spheres in the QM snapshots show uncertainty in the delocalized 
protons positions. For clarity all other particles are represented by their centroids. In the bottom panels, probabilities are 
shown on a log scale and are normalized by their maximum values. 


change in the ionization of Tyr57 obtained experimen¬ 
tally upon exchanging H for D (7%) is in good agreement 
with the value predicted from our simulations (4.1%). 
The under-prediction of the isotope effect on fractional 
ionization from our simulations is in line with the slightly 
low value of the simulated excess isotope effect, consis¬ 
tent with recent observations that the B3LYP-D3 density 
functional slightly underestimates the degree of proton 
sharing, and hence isotope effects, in hydrogen bonded 
systems upon including nuclear quantum effects 37 . 

The large degree of proton sharing with the depro- 
tonated Tyr57 residue upon including nuclear quantum 
effects can be elucidated by considering the potential en¬ 
ergy required, Ai% = o, to move a proton in the KSI ty¬ 
rosine triad from its energetic minimum to a perfectly 
shared position between the two tyrosine groups (y = 0). 
Ai% = o depends strongly on the positions of the residues 
comprising the triad, and in particular on the separation 
between the proton donor and acceptor oxygen atoms. 
Figure 4a shows Ai% = o computed as a function of the 
distance between 016 and 057, Roo , for the tyrosine 
triad in the absence of the protein environment. Re¬ 
moving the protein environment allows us to examine 
how changes in the triad distances from their positions 
in the enzyme affects the proton delocalization behavior 
without introducing steric overlaps with other active site 
residues (see SI Methods Id). 

Figure 4a shows that for the range of oxygen distances 
observed in the tyrosine triad ( Roo — 2.50 - 2.65 A) 


A£7=o for H16 is 3-6 kcal/mol. This is 6-12 times the 
thermal energy (fc^T) available at 300 K, leading to a 
very low thermal probability of the proton shared state 
(lower than e~ 6 = 2 x 10 -3 ). However, upon including 
nuclear quantum effects the system possesses ZPE, which 
in this system is ^4 kcal/mol. The ZPE closely matches 
AE U=0 and thus floods the potential energy wells along 
the proton sharing coordinate (Fig. 4b) allowing facile 
proton sharing (Fig. 3), i.e. inducing a transition to 
a LBHB type regime where the protons are quantum 
mechanically delocalized between the hydrogen bonded 
heavy atoms 18-20 . This leads to qualitatively different 
behavior of the protons in the active site of KSI D40Ar : 
from classical hydrogen bonding to quantum delocaliza¬ 
tion. The proton delocalization between the residues al¬ 
lows for ionization to be shared among three tyrosines 
to stabilize the deprotonation of Tyr57, leading to the 
large observed pK a shift relative to the value in the clas¬ 
sical limit (Fig. S2). This change in proton behavior 
gives rise to the large excess isotope effect since an O- 
D stretch possesses a ZPE of ~3 kcal/mol, which is no 
longer sufficient to fully flood the potential energy well in 
the proton sharing coordinate v. As the 0-0 separation 
is decreased below the values observed in KSI’s tyrosine 
triad, AF^=o becomes negligible compared to the ther¬ 
mal energy (^0.6 kcal/mol at 300 K). Hence, at very 
short distances (< 2.3 A) thermal fluctuations alone per¬ 
mit extensive proton sharing between the residues and 
the ZPE plays a negligible role in determining the pro- 
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FIG. 4. Comparison of the energy required to share a pro¬ 
ton between residues (AE u =o) as a function of the hydrogen 
bond donor-acceptor 0-0 distance ( Roo ) as compared to the 
zero-point energy, (a) AE u= o as a function of the oxygen- 
oxygen distance (Roo) between 016 and 057 using the ty¬ 
rosine triad geometry from a crystal structure (details in SI 
Methods Id). Inset figure shows the probability distribution 
of Roo obtained from the AI-PIMD simulation of KSI D40iV 
with ionized Tyr57. The probabilities are normalized by their 
maximum values, (b) Potential energy as a function of the 
proton transfer coordinate, u, for Roo = 2.6 A indicating 
values for the hydrogen and deuterium (0-H and 0-D) zero 
point energies, ZPEh and ZPEd • The position of Tyr32 is 
fixed as the proton H16 is scanned along Roo • 



FIG. 5. Simulation snapshot of the KSI D40iV active site with 
the bound intermediate analog (phenol) that gives rise to an 
extended delocalized proton network. Green, red and white 
represent C, O and H atoms. The blue spheres represent 
uncertainty in the delocalized proton position in the hydrogen 
bond network. For clarity all other particles are represented 
by their centroids. 


tons positions. Thus one would expect a small isotope 
effect. On the other hand, at bond lengths in excess 
of 2.7 A, A£^=o becomes so large (> 8 kcal/mol, Fig. 
4a) that the ZPE is not sufficient to flood the barriers, 
also resulting in a small expected isotope effect 40 . The 
large excess isotope effect in KSl D40N thus arises from 
the close matching of the ZPE and depth of the ener¬ 
getic well (AE u= q) which is highly sensitive to the 0-0 
distance. Hence although proton delocalization can oc¬ 
cur classically at short 0-0 distances (< 2.3 A), nuclear 
quantum effects allow this to occur for a much wider 
range of 0-0 distances (up to ^2.6 A), making delocal¬ 
ization feasible without incurring the steep steric costs 
that would be associated with bringing oxygen atoms any 
closer. The distances in KSI’s active site triad motif thus 
maximize quantum proton delocalization, which acts to 
stabilize the deprotonated residue. 


D. Quantum delocalization stabilizes an intermediate 
analog by sharing ionization 

Finally, we considered the role of nuclear quantum ef¬ 
fects when an intermediate analog participates in the 
active site hydrogen bond network. Recent experi¬ 
ments have investigated how the binding of intermediate 
analogs to KSI D40Ar affects the sharing of ionization along 
the extended hydrogen bond network that is formed (Fig. 
lb) 15 . These experiments identified that ionization shar¬ 
ing is maximized when phenol, whose solution pK a of 10 
equals that of the actual intermediate of KSI (Fig. la), is 
bound 15,41,42 . We thus performed AI-PIMD simulations 


of the KSI D40Ar -phenol complex. The protons were ob¬ 
served to be delocalized across the network (Fig. 5) with 
partial ionizations of Tyr57, Tyrl6 and phenol calculated 
to be 18.5%, 56.7% and 22.5% (± 0.7%), compared with 
estimates from previous experiments using 13 C NMR of 
40%, 40%, 20% 15 . Hence simulation and experiment are 
in good agreement that the ionization is shared almost 
equally among the three residues, i.e. that there is al¬ 
most no difference (< 2/c^T) in the free energy upon 
shifting the ionization among any of the three groups. 
Therefore, the ability of protons to delocalize within the 
KSi^ 40Ar tyrosine triad, initially found in the enzyme in 
the absence of the intermediate analog and manifested 
as a strongly perturbed pK a , extends upon incorpora¬ 
tion of the intermediate analog, which shifts the center 
of the ionization along the network from Tyr57 to Tyrl6. 
In both cases, proton delocalization acts to share the ion¬ 
ization of a negatively charged group, which suggests that 
KSI could utilize quantum delocalization in its active site 
hydrogen bond network to distribute the ionization aris¬ 
ing in its intermediate complex (Fig. la) so as to provide 
energetic stabilization. 


III. CONCLUSIONS 

In conclusion, KSI D40Ar exhibits a large equilibrium 
isotope effect in the acidity of its active site tyrosine 
residues arising from a highly specialized triad motif, con¬ 
sisting of several short 0-0 distances, whose positions 
enhance quantum delocalization of protons within the 
active site hydrogen bond network. This delocalization 
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manifests in a very large isotope effect and substantial 
acidity shift. Our simulations, which include electronic 
quantum effects and exactly treat the quantum nature of 
the nuclei, show qualitatively and quantitatively different 
proton behavior compared to conventional simulations in 
which the nuclei are treated classically, and provide good 
agreement with experiment. The ability to perform such 
simulations thus offers the opportunity to investigate in 
unprecedented detail the plethora of systems in which 
short-strong hydrogen bonds occur, where incorporating 
both nuclear and electronic quantum effects is crucial to 
understand their functions. 


IV. MATERIALS AND METHODS 

A. Expression and purification of KSI 

Wild-type and D40N KSI from Pseudomonas putida 
were over-expressed in BL-21 A1 cells (Invitrogen), 
isolated by affinity chromatography using a custom- 
designed deoxycholate-bound column resin, and purified 
by gel-filtration chromatography (GE Healthcare) as de¬ 
scribed previously 43 . For 13 C NMR experiments, 13 C^- 
tyrosine was incorporated into KSI according to methods 
described previously 15 . 


B. UV-Vis titration experiments 

A series of buffers was prepared with a pL between 
4 and 10 by weighing portions of a weak acid and its 
sodium-conjugate base salt, and adding the appropriate 
form of distilled deionized water (Millipore H 2 O, Spectra 
stable isotopes sterile-filtered D 2 0 (>99% 2 H)). Buffers 
were prepared at 40 mM. Tyr57 is solvent accessible, so 
the tyrosine residues in the active site network are ex¬ 
pected to be fully deuterated in D 2 0 solution. 

The following buffer systems were used for the fol¬ 
lowing pL ranges: acetic acid/sodium acetate, 4-5.25; 
sodium monobasic phosphate/dibasic phosphate, 5.5- 
8.25; sodium bicarbonate/sodium carbonate, 8.5-10. 
Buffers were stored at room temperature with caps firmly 
sealed. 

After preparation of buffers, pL was recorded using 
an Orion2 Star glass electrode (Thermo), immediately 
following calibration with standard buffers at pH 4, 7, 
and 10. In H 2 0, the pH of the buffer was taken as the 
reading on the electrode. In D 2 0, the pD of the buffer 
was calculated by adding 0.41 to the operational pH* 
from the electrode reading 44 . A series of samples for 
titration was prepared by combining 60 fiL protein (100 
//M stock in buffer-free L 2 0), buffer (150 /iL of 40 mM 
stock), and extra L 2 0. The final samples were 600 /iL, 
10 /iM protein, 10 mM buffer. 

UV-Vis measurements were carried out on the sam¬ 
ples on a Lambda 25 spectrophotometer (Perkin Elmer), 


acquiring data from 400 to 200 nm with a 1.0 nm data in¬ 
terval and 960 nm/min scan rate and 1.00 nm slit width. 
For each measurement a background was taken to the 
pure buffer of a given pL, before acquiring on the protein- 
containing sample. Spectra were recorded in duplicate to 
control for random detector error. 

The spectra were baselined by setting the absorption at 
320 nm to zero, and the change in absorption at 300 nm 
was followed at varying pLs, employing a previous estab¬ 
lished method 16 to determine the fractional ionization of 
a tyrosine-tyrosinate pair. The error in A 300 from com¬ 
paring duplicate spectra following baselining was gen¬ 
erally between 0-2%. For each pL, the average A 300 
was calculated and converted to an extinction coefficient 
( 6300 )- The titration experiment was repeated on two in¬ 
dependently prepared buffer stocks to control for error in 
buffer preparation. 


C. Simulations 

AI-PIMD and AIMD simulations were performed using 
a QM/MM approach of KSI D40iV with Tyr57 protonated, 
KSi^ 40Ar with Tyr57 ionized, KSI 1140 ^ with the inter¬ 
mediate analog bound, and tyrosine in aqueous solution. 
The simulations were carried out in the NVT ensemble 
at 300 K with a time step of 0.5 fs. The path integral 
generalized Langevin equation (PIGLET) approach was 
utilized, which allowed results within the statistical error 
bars to be obtained using only 6 path integral beads to 
represent each particle 34 . The electronic structure in the 
QM region was evaluated using the B3LYP functional 28 
with dispersion corrections 29 . The 6-31G* basis set was 
used as we found it to produce proton transfer poten¬ 
tial energy profiles with a mean absolute error of less 
than 0.4 kcal/mol for this system (see Fig. S6). Energies 
and forces in the QM region and the electrostatic interac¬ 
tions between the QM and MM regions were obtained us¬ 
ing an MPI interface to the GPU-accelerated TeraChem 
package 35,36 . Atoms in the MM region were described 
using the AMBER03 force field 45 , and the TIP3P wa¬ 
ter model 46 . The simulations were performed using peri¬ 
odic boundary conditions with Ewald summation to treat 
long-range electrostatic interactions. The energies and 
forces within the MM region, and the Lennard-Jones in¬ 
teractions between the QM and MM regions were calcu¬ 
lated by MPI calls to the LAMMPS molecular dynamics 
package 47 . The QM region of KSI D40Ar contained the p- 
methylene phenol side chains of residues Tyrl6, Tyr32 
and Tyr57 (Fig. S5a). For KSL 040 ^ with the intermedi¬ 
ate analog, residue Aspl03 and the bound intermediate 
analog were also included in the QM region (Fig. S5b). 
The QM region of tyrosine in solution contained the side 
chain of the tyrosine residue and the 41 water molecules 
within 6.5 A of the side chain O-H group. All bonds 
across the QM/MM interface were capped with hydro¬ 
gen link atoms in the QM region 26 . These capping atoms 
were constrained to be along the bisected bonds and do 
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not interact with the MM region. 

The initial configuration of KSI D40Ar was obtained 
from a crystal structure 5 (PDB ID 10GX). For KSI D40iV 
with the intermediate analog, a crystal structure 15 (PDB 
ID 3VGN) was used with the ligand changed to phe¬ 
nol. The crystal structures were solvated in TIP3P water 
and energy minimized before performing AI-PIMD sim¬ 
ulations. The initial configuration for tyrosine in aque¬ 
ous solution was obtained by solvating the amino acid 
in TIP3P 46 water using the AMBER03 force field 45 and 
equilibrating for 5 ns in the NPT ensemble at a temper¬ 
ature of 300 K and pressure of 1 bar. Each system was 
then equilibrated for 10 ps followed by production runs 
of 30 ps. 

To calculate the excess isotope effect, A A pK a , (SI 
Methods Ic) we used the thermodynamic free energy per¬ 
turbation (TD-FEP) path integral estimator 22 . Com¬ 
bined with an appropriate choice of the integration vari¬ 
able to smooth the free energy derivatives 22 , this allowed 
us to evaluate the isotope effects in the liquid phase us¬ 
ing only a single AI-PIMD trajectory. Simulations per¬ 
formed with D substitution showed no change within the 
statistical error bars reported. 
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